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Abstract 

Ultrasound computed tomography (USCT) holds great promise for improving the detection and 
management of breast cancer. Because they are based on the acoustic wave equation, waveform inversion- 
based reconstruction methods can produce images that possess improved spatial resolution properties 
over those produced by ray-based methods. However, waveform inversion methods are computationally 
demanding and have not been applied widely in USCT breast imaging. In this work, source encoding 
concepts are employed to develop an accelerated USCT reconstruction method that circumvents the 
large computational burden of conventional waveform inversion methods. This method, referred to as 
the waveform inversion with source encoding (WISE) method, encodes the measurement data using 
a random encoding vector and determines an estimate of the sound speed distribution by solving a 
stochastic optimization problem by use of a stochastic gradient descent algorithm. Both computer- 
simulation and experimental phantom studies are conducted to demonstrate the use of the WISE method. 

The results suggest that the WISE method maintains the high spatial resolution of waveform inversion 
methods while significantly reducing the computational burden. 
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I. Introduction 

After decades of research m-gi, advancements in hardware and computing technologies are 
now facilitating the clinical translation of ultrasound computed tomography (USCT) for breast 
imaging applications |21, GO-GOl. USCT holds great potential for improving the detection and 
management of breast cancer since it provides novel acoustic tissue contrasts, is radiation- and 
breast-compression-free, and is relatively inexpensive. 0, lITOll . Several studies have reported 
the feasibility of USCT for characterizing breast tissues El, 10-I61, ifTOl . iflTl . Although some 
USCT systems are capable of generating three images that depict the breast’s acoustic reflectivity, 
acoustic attenuation, and sound speed distributions, this study will focus on the reconstruction 
of the sound speed distribution. 

A variety of USCT imaging systems have been developed for breast sound speed imaging 
lf5l , 1171 , m, lfT2l - lfl5l . In a typical USCT experiment, acoustic pulses that are generated by 
different transducers are employed, in turn, to insonify the breast. The resulting wavefield data 
are measured by an array of ultrasonic transducers that are located outside of the breast. Here 
and throughout the manuscript, a transducer that produces an acoustic pulse will be referred 
to as an emitter; the transducers that receive the resulting wavefield data will be referred to as 
receivers. From the collection of recorded wavefield data, an image reconstruction method is 
utilized to estimate the sound speed distribution within the breast 0, 0, lUOil . 

The majority of USCT image reconstruction methods for breast imaging investigated to date 
have been based on approximations to the acoustic wave equation lH2l . [il6l - [l24l . A relatively 
popular class of methods is based on geometrical acoustics, and are commonly referred to as ‘ray- 
based’ methods. These methods involve two steps. First, time-of-flight (TOF) data corresponding 
to each emitter-receiver pair are estimated E51 . Under a geometrical acoustics approximation, 
the TOF data are related to the sound speed distribution via an integral geometry, or ray-based, 
imaging model Ifl6l . Il26l Second, by use of the measured TOF data and the ray-based imaging 
model, a reconstruction algorithm is employed to estimate the sound speed distribution. Although 
ray-based methods can be computationally efficient, the spatial resolution of the images they 
produce is limited due to the fact that diffraction effects are not modeled |[23l [[271 This is 
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undesirable for breast imaging applications, in which the ability to resolve fine features, e.g., 
tumor spiculations, is important for distinguishing healthy from diseased tissues. 

USCT reconstruction methods based on the acoustic wave equation, also known as full-wave 
inverse scattering or waveform inversion methods, have also been explored for a variety of 
applications including medical imaging llT2l . l(22l . Ii23ll . Il28ll and geophysics |[29il - lf3Tfl . Because 
they account for higher-order diffraction effects, waveform inversion methods can produce images 
that possess higher spatial resolution than those produced by ray-based methods ll23l . Il28l . 
However, conventional waveform inversion methods are iterative in nature and require the wave 
equation to be solved numerically a large number of times at each iteration. Consequently, 
such methods can be extremely computationally burdensome. For special geometries llT2l . Il32ll . 
efficient numerical wave equation solvers have been reported. However, apart from special cases, 
the large computational burden of waveform inversion methods has hindered their widespread 
application. 

A natural way to reduce the computational complexity of the reconstruction problem is to 
reformulate it in a way that permits a reduction in the number of times the wave equation 
needs to be solved. In the geophysics literature, source encoding methods have been proposed 
to achieve this 11291 - 1(311 . When source encoding is employed, at each iteration of a prescribed 
reconstruction algorithm, all of the acoustic pulses produced by the emitters are combined (or 
‘encoded’) by use of a random encoding vector. The measured wavefield data are combined in 
the same way. As a result, the wave equation may need to be solved as few as twice at each 
algorithm iteration. In conventional waveform inversion methods, this number would be equal 
to twice the number of emitters employed. Although conventional waveform inversion methods 
may require fewer algorithm iterations to obtain a specified image accuracy compared to source 
encoded methods, as demonstrated later, the latter can greatly reduce the overall number of times 
the wave equation needs to be solved. 

In this study, a waveform inversion with source encoding (WISE) method for USCT sound 
speed reconstruction is developed and investigated for breast imaging with a circular trans¬ 
ducer array. The WISE method determines an estimate of the sound speed distribution by 
solving a stochastic optimization problem by use of a stochastic gradient descent algorithm Il30ft . 
[f33ll . Unlike previously studied waveform inversion methods that were based on the Helmholtz 
equation 1(22(1 . l(23l . the WISE method is formulated by use of the time-domain acoustic wave 
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equation [f34l - [|36ll and utilizes broad-band measurements. The wave equation is solved by use 
of a computationally efficient k-space method that is accelerated by use of graphics processing 
units (GPUs). In order to mitigate the interference of the emitter on its neighboring receivers, a 
heuristic data replacement strategy is proposed. The method is validated in computer-simulation 
studies that include modeling errors and other physical factors. The practical applicability of the 
method is further demonstrated in studies involving experimental breast phantom data. 

The remainder of the paper is organized as follows. In Section II, USCT imaging models in 
their continuous and discrete forms are reviewed. A conventional waveform inversion method and 
the WISE method for sound speed reconstruction are formulated in Section III. The computer- 
simulation studies and corresponding numerical results are presented in Sections IV and V, 
respectively. In Section VI, the WISE method is further validated in experimental breast phantom 
studies. Finally, the paper concludes with a discussion in Section VII. 

II. Background: USCT imaging models 

In this section, imaging models that provide the basis for image reconstruction in waveform 
inversion-based USCT are reviewed. 


A. USCT imaging model in its continuous form 


Although a digital imaging system is properly described as a continuous-to-discrete (C-D) 
mapping (See Chapter 7 in If37l ). for simplicity, a USCT imaging system is initially described 
in its continuous form below. 

In USCT breast imaging, a sequence of acoustic pulses is transmitted through the breast. 
We denote each acoustic pulse by s m (r,t) £ Lr(M 3 x [0, oo)), where each pulse is indexed by 
an integer m for m = 0,1, • • • , M — 1 with M denoting the total number of acoustic pulses. 
Although it is spatially localized at the emitter location, each acoustic pulse can be expressed as 
a function of space and time. When the m-th pulse propagates through the breast, it generates a 
pressure wavefield distribution denoted by p m ( r, t) e L 2 (R 3 x [0, oo)). If acoustic absorption and 
mass density variations are negligible, p m (r,t) in an unbounded medium satisfies the acoustic 
wave equation lf38l : 


V 2 p m (r,t) 


1 d 2 

c 2 (r) dt 2 


Pm(r,t) 


-47rs m (r,t), 


( 1 ) 
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where c(r) is the sought-after sound speed distribution. Equation © can be expressed in operator 
form as 

Pm{r,t) =n c s m (r,t), (2) 

where the linear operator TL C : L 2 (R 3 x [0, oo)) i— > L 2 (R 3 x [0, oo)) denotes the action of the wave 
equation and is independent of the index of m. The superscript ‘c’ indicates the dependence of 
TL C on c(r). 

Consider that p m (r,t) is recorded outside of the object for r G Vt m and t G [0, T], where 
Q m c R 3 denotes a continuous measurement aperture. In this case, when discrete sampling 
effects are neglected, the imaging model can be described as a continuous-to-continuous (C-C) 
mapping as: 

9 m(r,t) = M m n c s m (r,t), for m = 0,1, • ■ ■ , M - 1, (3) 

where g m (r,t) G L 2 (Q rn x [0, T]) denotes the measured data function and the operator A4 m is 
the restriction of TL C to Q m x [0, T], The m-dependent operator A4 rn allows Eqn. © to describe 
USCT imaging systems in which the measurement aperture varies with emitter location. Here 
and throughout the manuscript, we will refer to the process of firing one acoustic pulse and 
acquiring the corresponding wavefield data as one data acquisition indexed by m. The USCT 
reconstruction problem in its continuous form is to estimate the sound speed distribution c(r) 
by use of Eqn. © and the data functions {g m (r,t)}^lQ. 

B. USCT imaging model in its discrete forms 

A digital imaging system is accurately described by a continuous-to-discrete (C-D) imaging 
model, which is typically approximated in practice by a discrete-to-discrete (D-D) imaging model 
to facilitate the application of iterative image reconstruction algorithms. A C-D description of the 
USCT imaging system is provided in Appendix lAl Below, a D-D imaging model for waveform- 
based USCT is presented. This imaging model will be employed subsequently in the development 
of the WISE method in Section |]T[1 

Construction of a D-D imaging model requires the introduction of a finite-dimensional approx¬ 
imate representations of the functions c(r) and s m (r, i), which will be denoted by the vectors 
c G R ;Y and s m G R AV '. Here, N and L denote the number of spatial and temporal samples, 
respectively, employed by the numerical wave equation solver. In waveform-based USCT, the 
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way in which c(r) and s m (r, t) are discretized to form c and s m is dictated by the numerical 
method employed to solve the acoustic wave equation. In this study, we employ a pseudospectral 
k-space method Il34ll - ll36ll . Accordingly, c(r) and s m (r,t) are sampled on Cartesian grid points 
as 

[c] n = c(r n ), and [s m ] nL+l = s m (r n , /A*), for ( 4 ) 

where A* denotes the temporal sampling interval and r n denotes the location of the n-th point. 

For a given c and s m , the pseudospectral k-space method can be described in operator form 
as 

p a m = H c s m , (5) 

where the matrix H c is of dimension NL x NL and represents a discrete approximation of the 
wave operator LL C defined in Eqn. ©, and the vector p', l „ represents the estimated pressure data 
at the grid point locations and has the same dimension as s m . The superscript ‘a’ indicates that 
these values are approximate, i.e., \p' rr ] n i.+i ~ Pm.(j n , /A 4 ). We refer the readers to If34l - ll36l for 
additional details regarding the pseudospectral k-space method. 

Because the pseudospectral k-space method yields sampled values of the pressure data on a 
Cartesian grid, a sampling matrix M m is introduced to model the USCT data acquisition process 
as 

g a m = M m p* m = M m H c s m . (6) 


Here, the N rec L x NL sampling matrix M m extracts the pressure data corresponding to the 
receiver locations on the measurement aperture Q m , with N rec denoting the number of receivers. 
The vector g a t denotes the predicted data that approximates the true measurements. In principle, 
M m can be constructed to incorporate transducer characteristics, such as finite aperture size and 
temporal delays. For simplicity, we assume that the transducers are point-like in this study. When 
the receiver and grid point locations do not coincide, an interpolation method is required. As an 
example, when a nearest-neighbor interpolation method is employed, the elements of M m are 
defined as 


[M m ] 


n Tec L+l,nL+l 


1, for n = l m (n Tec ), 
0, otherwise, 


(7) 


where [’ ^- m \n Tec L+i,nL+i denotes the element of M m at the (n rec L + f)-th row and the (nL + l )-th 
column, and l rn (n rec ) denotes the index of the grid point that is closest to r(m, n rec ). Here, 
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r (m, n rec ) denotes the location of the n rec -th receiver in the m-th data acquisition. Equation ([6]) 
represents the D-D imaging model that will be employed in the remainder of this study. 


III. Waveform inversion with source encoding for USCT 


A. Sequential waveform inversion in its discrete form 


A conventional waveform inversion method that does not utilize source encoding will be 
employed as a reference for the developed WISE method and is briefly described below. Like 
other conventional approaches, this method sequentially processes the data acquisitions g m for 
m — 0,1, ■ • • , M — 1 at each iteration of the associated algorithm. As such, we will refer to the 
conventional method as a sequential waveform inversion method. 

A sequential waveform inversion method can be formulated as a non-linear numerical opti¬ 
mization problem: 

c = argmin{J r (c) + /371(c)}, (8) 

C 


where ^(c), IZ(c), and 3 denote the data fidelity term, the penalty term, and the regularization 
parameter, respectively. The data fidelity term J-(c) is defined as a sum of squared £ 2 -norms of 
the data residuals corresponding to all data acquisitions as: 

1 M—l 

T(c) = - llgm - M m H C S m || 2 , (9) 

m =0 

where g m 6 L denotes the measured data vector at the m-th data acquisition. The choice 
of the penalty term will be addressed in Section [TV] 

The gradient of J r (c) with respect to c, denoted by J, will be computed by discretizing an 
expression for the Frechet derivative that is derived assuming a continuous form of Eqn. ©. 
The Frechet derivative is described in Appendix [B] Namely, the gradient is approximated as 


[JL = 


M—l 

El 

m =0 


M—l L—2 


wrEE 

L In m= o i = i 


QmJ nL-\-(L—l) 


feU+Fl - 2[p UnL+l + [Pm]»I+i+l 
A* 


( 10 ) 


where J m denotes the gradient of |||g m — M m H c s m || 2 with respect to c and the vector 
contains samples that approximate adjoint wavefield q m (r, t) that satisfies Eqn. (l34l) in Appendix 
|Bl By use of the pseudospectral k-space method, q),, can be calculated as 


q a 



mi 


(ID 
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Here, N m 


[" 7 ” m\nL+l 

{n : X m (n rec ),n rec = 


[gm - &»W(n)i+M’ if n e N m , 


( 12 ) 


0, otherwise 

0,1, ■ ■ ■ , N rec — 1}, and I ~ 1 denotes the inverse mapping of 


- L 'm • 

Given the explicit form of J in Eqn. (fTOl) . a variety of optimization algorithms can be employed 
to solve Eqn. ([8]) [[39]]. Algorithm 1 describes a gradient descent-based sequential waveform 
inversion method. Note that at every algorithmic iteration, the sequential waveform inversion 
method updates the sound speed estimate only once using the gradient J accumulated over all 
J m for m — 0,1, • • • , M — 1. This is unlike the Kaczmarz method—also known as the algebraic 
reconstruction technique lU6ll . fl9| . fl40ll —that updates the sound speed estimate multiple times 
in one algorithmic iteration. In Line-10 of Algorithm 1, J R denotes the gradient of 71(c) with 
respect to c. 


Algorithm 1 Gradient descent-based sequential waveform inversion. 
Input: {g™,}, {s m }, cA 0) 

Output: c 

l: k <— 0 {k is the number of algorithm iteration.} 

2 : while stopping criterion is not satisfied do 

3: k <— k + 1 

4: J <— 0 

5: for m := 0 to M — 1 do 

6: •*— H c s m {m is the index of the emitter.} 

7: q} n -t— H c r m (r m is calculated via Eqn. (fill).} 

8: J J + J m {J m is calculated via Eqn. (fTOl).} 

9: end for 

10: J <- J + /3J R 

li: Determine step size A via a line search 

12: C •<— C (fc-1) — AJ 

13: end while 

14: C = 
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In Algorithm [H H c is the most computationally burdensome operator, representing one run 
of the wave equation solver. Note that it appears in Lines-6, -7, and -11. Because Lines-6 and 
-7 have to be executed M times to process all of the data acquisitions, the wave equation solver 
has to be executed at least (2M + 1) times at each algorithm iteration. The line search in Line-11 
searches for a step size along the direction of — J so that the cost function is reduced by use of 
a classic trial-and-error approach ll39l . Note that, in general, the line search will require more 
than one application of H c , so (2 M + 1) represents a lower bound on the total number of wave 
equation solver runs per iteration. 

B. Stochastic optimization-based waveform inversion with source encoding (WISE) 

In order to alleviate the large computational burden presented by sequential waveform inversion 
methods (e.g., Algorithm 1), a source encoding method has been proposed ll22l . ||29l , lf41~ll . This 
method has been formulated as a stochastic optimization problem and solved by various stochastic 
gradient-based algorithms l(30l . EHl . In this section, we adapt the stochastic optimization-based 
formulation in ll30l to find the solution of Eqn. d8]). 


Algorithm 2 Waveform inversion with source encoding (WISE) algorithm. 

Input: {gm}, {s m }, c (0) 

Output: c 

l: k <— 0 {k is the number of algorithm iteration} 

2 : while stopping criterion is not satisfied do 

3: k k + 1 

4: Draw elements of w from independent and identical Rademacher distribution. 

5: p w H c s w {s w is calculated via Eqn. (fl4l) . } 

6: q w H't w { t w is calculated via Eqn. (fT71).} 

7: J •<— J w + /3J R {J w is calculated via Eqn. (fT6l)} 

8: Determine step size A by use of line search 

9: C <— C (fc-1) — AJ 

to: end while 

11: C = C 
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The WISE method seeks to minimize the same cost function as the sequential waveform 
inversion method, namely, Eqn. ([8]). However, to accomplish this, the data fidelity term in Eqn. ([9]) 
is reformulated as the expectation of a random quantity as 112911 - 113111 . (331 , (4T1 . Ii42il 

T s (c) = E w {l||g w - MH c s w || 2 }, (13) 

where E w denotes the expectation operator with respect to the random source encoding vector 
w e R M , M = M m is the sampling matrix that is assumed to be identical for m = 0,1, • • • , M— 
1, and g w and s w denote the w-encoded data and source vectors, defined as 

M—l M—1 

g w = y>] ro6m , and s w = y^[w] m s m , (14) 

m= 0 m=0 

respectively. It has been demonstrated that Eqns. © and (Il3l) are mathematically equivalent 
when w possesses a zero mean and an identity covariance matrix (30l . (33l . (421. In this case, 
the optimization problem whose solution specifies the sound speed estimate can be re-expressed 
in a stochastic framework as 

c = argmin E w {A||g w - MH c s w || 2 } + pn(c), (15) 

c ^2 


which we refer to as the waveform inversion with source encoding (WISE) method. An im¬ 
plementation of the WISE method that utilizes the stochastic gradient descent algorithm is 
summarized in Algorithm [2} 

In Algorithm [2] the wave equation solver needs to be run one time in each of Lines-5 and 
6. In the line search to determine the step size in Line 8, the wave equation solver needs to be 
run at least one time, but in general will require a small number of additional runs, just as in 
Algorithm Q] Accordingly, the lower bound on the number of required wave equation solver runs 
per iteration is 3, as opposed to (2 M + 1) for the conventional sequential waveform inversion 
method described by Algorithm 1. As demonstrated in geophysics applications (291 . (311 . (4T1 
and the breast imaging studies below, the WISE method provides a substantial reduction in 
reconstruction times over use of the standard sequential waveform inversion method. In Line-7, 
J w can be calculated analogously to Eqn. (UOl) as 


L—2 

E 

i=i 


Q W ]nL+(L-l) 


[p W ]nL+Z—1 


2 [p W ]nL+l + [p W ]nL+Z+l 


A* 


(16) 
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Here, we drop the subscript m of both X x (n) and N because we assume M to be identical for 



(17) 


all data acquisitions. Various probability density functions have been proposed to describe the 
source encoding vector w |[29l . ETA . BT1 . In this study, we employed a Rademacher distribution 
as suggested by ll29l . in which case each element of w had a 50% chance of being either +1 
or —1. 


IV. Description of computer-simulation studies 


Two-dimensional computer-simulation studies were conducted to validate the WISE method 
for breast sound speed imaging and demonstrate its computational advantage over the standard 
sequential waveform inversion method. 

A. Measurement geometry 

A circular measurement geometry was chosen to emulate a previously reported USCT breast 
imaging system iflOl . If23l . Il43l . As depicted in Fig. [H 256 ultrasonic transducers were uniformly 
distributed on a ring of radius 110 mm. The generation of one USCT data set consisted of M = 
256 sequential data acquisitions. In each data acquisition, one emitter produced an acoustic pulse. 
The acoustic pulse was numerically propagated through the breast phantom and the resulting 
wavefield data were recorded by all transducers in the array as described below. Note that the 
location of the emitter in every data acquisition was different from those in other acquisitions, 
while the locations of receivers were identical for all acquisitions. 

B. Numerical breast phantom 

A numerical breast phantom of diameter 98 mm was employed. The phantom was composed 
of 8 structures representing adipose tissues, parenchymal breast tissues, cysts, benign tumors, 
and malignant tumors, as shown in Fig. [3 For simplicity, the acoustic attenuation of all tissues 
was described by a power law with a fixed exponent y — 1.5 ll44ll . The corresponding sound 
speed and the attenuation slope values are listed in TABFE U [|44il - [|46l . Both the sound speed 
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and the attenuation slope distributions in Fig. [2] were sampled on a uniform Cartesian grid with 
spacing A s = 0.25 mm. The finest structure (indexed by 7 in Fig. ERa)) was of diameter 3.75 
mm. 

C. Simulation of the measurement data 

1) First-order numerical wave equation solver: Acoustic wave propagation in acoustically 
absorbing media was modeled by three coupled first-order partial differential equations |[47l : 

J^ u (r,f) = —Vp(r,f) (18a) 

— p(r, t) = —V • u(r, t) + 47t / dt's(r,t') (18b) 

ut J o 

p(r,t) = c 2 (r) [l + r(r)-^-(—V 2 ) y/2 ~ 1 + 77(r)(-V 2 ) (y+1)/2_1 ]p(r, t), (18c) 

where u(r,t), p(r,t), and p(r) denote the acoustic particle velocity, the acoustic pressure, and 
the acoustic density, respectively. The functions r(r) and r/(r) describe acoustic absorption and 
dispersion during the wave propagation 11471 : 

r(r) = -2a 0 (r)c 0 (r) y_1 , rj( r) = 2a 0 (r)c 0 (r) y tan(vn//2), (19) 

where a 0 (r) and y are the attenuation slope and the power law exponent, respectively. When the 
medium is assumed to be lossless, i.e., a 0 (r) — 0, it can be shown that Eqn. (1181) is equivalent 
to Eqn. (IB- 

Based on Eqn. (1181) . a pseudospectral k-space method was employed to simulate acoustic 
pressure data ll36l . 11471 . This method was implemented by use of a first-order numerical scheme 
on GPU hardware. The calculation domain was of size 512 x 512 mm 2 , sampled on a 2048 x 2048 
uniform Cartesian grid of spacing A s = 0.25 mm. A nearest-neighbor interpolation was employed 
to place all transducers on the grid points. On a platform consisting of dual quad-core CPUs with 
a 3.30 GHz clock speed, 64 gigabytes (GB) of random-accessing memory (RAM), and a single 
NVIDIA Tesla K20 GPU, the first-order pseudospectral k-space method required approximately 
108 seconds to complete one forward simulation. 

2) Acoustic excitation pulse: The excitation pulse employed in this study was assumed to be 
spatially localized at the emitter location while temporally it was a f c — 0.8 MHz sinusoidal 
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function tapered by a Gaussian kernel with standard deviation a = 0.5 ps, i.e., 

{exp ( — ) sin(27r f c t), at the m-th emitter location 

' w ^ (20) 

0, otherwise, 

where the constant time shift t c = 3.2 /is. The temporal profile and the amplitude frequency 
spectrum of the excitation pulse are plotted in Fig. [3]-(a) and -(b), respectively. The excitation 
pulse contained approximately 3 cycles. 

3) Generation of non-attenuated and attenuated noise-free data: For every data acquisition 
(indexed by m), the first-order pseudospectral k-space method was run for 3600 time steps with 
a time interval A* = 0.05 /is (corresponding to a 20 MHz sampling rate). Downsampling the 
recorded data by taking every other time sample resulted in a data vector g m (see Eqn. (HJ) that 
was effectively sampled at 10 MHz and was of dimensions ML with M = 256 and L = 1800. 
The data vector at the 0-th data acquisition, g 0 , is displayed as a 2D image in Fig. IU-(a). This 
undersampling procedure was introduced to avoid inverse crime l|48ll so that the data generation 
and the image reconstruction employed different numerical discretization schemes. Repeating the 
calculation for m — 0,1, • • ■ , 255, we obtained a collection {g m } of data vectors that together 
represented one complete data set. Utilizing the absorption phantom described in Section IIV-B1 
a complete attenuated data set was computed. An idealized, non-attenuated, data set was also 
computed by setting a 0 (r) = 0. 

4) Generation of incomplete data: An incomplete data set in this study corresponds to one 
in which only N Tec receivers located on the opposite side of the emitter record the pressure 
wavefield, with N rec < M. Taking the 0-th data acquisition as an example (see Fig. |T|), only 
N rec = 100 receivers, indexed from 78 to 177, record the wavefield, while other receivers record 
either unreliable or no measurements. Incomplete data sets formed in this way can emulate two 
practical scenarios: (1) Signals recorded by receivers near the emitter are unreliable and therefore 
discarded lf23ll : and (2) An arc-shaped transducer array is employed that rotates with the emitter 

m, m, m. 

Specifically, incomplete data sets were generated as 


[g“ pl ] 


n Tec L+l 


[Sm] 


J m (n™ c )L+V 


for m =0,l. - ,M—1 
IUI n r ec = 01 i ... )jv rec_ 1 , 


( 21 ) 


where g™ cpl is the incomplete m-th data acquisition, which is of dimensions N rec L, with N rec < 
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M. The index map J m : {0,1, • • • , N rec — 1} i—^ M^ od is defined as 

J m (n rec ) = (m + n rec +---J mod M, (22) 

where ( m' mod M ) calculates the remainder of vn! divided by M, and the index set M^ >od 
collects indices of transducers that reliably record data at the m-th data acquisition and is defined 
as 

M^ od = | k mod M\k G [m + (M - N iec )/2 , m + (M + iV rec )/2) (23) 

Here, for simplicity, we assume that M and N rec are even numbers. In this study, we empirically 
set N Tec = 100 so that the object can be fully covered by the fan region as shown in Fig. Q] 

5) Generation of noisy data: An additive Gaussian white noise model was employed to 
simulate electronic measurement noise as 

gm = gm + n, (24) 

where g m and n are the noisy data vector and the Gaussian white noise vector, respectively. In 
this study, the maximum value of the pressure received by the 128-th transducer at the 0-th data 
acquisition with a homogeneous medium (water tank) was chosen as a reference signal amplitude. 
The noise standard deviation was set to be 5% of this value. An example of a simulated noiseless 
and noisy data acquisition is shown Fig. |U 

D. Image reconstruction 

1) Second-order pseudospectral k-space method: In the reconstruction methods described 
below, the action of the operator H c (Eqn. ([5])) was computed by solving Eqn. CO) by use of a 
second-order pseudospectral k-space method. This was implemented using GPUs. The calculation 
domain was of size 512 x 512 mm 2 , sampled on a 1024 x 1024 uniform Cartesian grid of spacing 
A s = 0.5 mm for reconstruction. On a platform consisting of dual octa-core CPUs with a 2.00 
GHz clock speed, 125 GB RAM, and a single NVIDIA Tesla K20C GPU, the second-order 
k-space method required approximately 7 seconds to complete one forward simulation. 

2) Sequential waveform inversion: To serve as a reference for the WISE method, we imple¬ 
mented the sequential waveform inversion method described in Algorithm 1. No penalty term 
was included (/3 = 0) because, due to its extreme computational burden, we only investigated this 
method in preliminary studies involving noise-free non-attenuated data. A uniform sound speed 
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distribution was employed as the initial guess, which corresponded to the known background 
value of 1.5 mm///s. The object was contained in a square region-of-interest (ROI) of dimension 
128 x 128 mm 2 (See Fig. [D, which was covered by 256 x 256 pixels. 

3) WISE method: We implemented the WISE method by use of Algorithm 2. Two types of 
penalties were employed in this study: a quadratic penalty expressed as 


^ Q ( c ) = 5 ^([ c W-h “ t c W+i-i) 2 + ([c]jjv«+i - [ c ] (j-i)N x +i) 2 , (25) 

3 * 

where N x and N y denote the number of grid points along the ‘x’ and ‘y’ directions respectively, 
and a total variation (TV) penalty, defined as |[50il . [15111 


^ TV ( C ) — ^ ^ \/ e + (MlVc+i — [ C ]jN x +i- 1) 2 + ([c]jN x +i — [c] , (26) 

3 * 

where e is a small number introduced to avoid dividing by 0 in the gradient calculation. In 
this study, we empirically selected e = 10 -8 . This value was fixed because we observed that 
it had a minor impact on the reconstructed images compared to the impact of f3. The use of 
this parameter can be avoided when advanced optimization algorithms are employed ll52ll . Il53ll . 
As in the sequential waveform inversion case, it was assumed that the background sound speed 
was known and the object was contained in a square ROI of dimension 128 x 128 mm 2 (See 
Fig. [Q). which corresponded to 256 x 256 pixels. The regularization parameters corresponding to 
the quadratic penalty and the TV penalty will be denoted by {3® and /3 TV , respectively. Optimal 
regularization parameter values should ultimately be identified by use of task-based measures of 
image quality ll37l . In this preliminary study, we investigated the impact of /! c - and /3 TV on the 
reconstructed images by sweeping their values over a wide range. 

4) Reconstruction from incomplete data: Because the WISE method requires M m to be 
identical for all m' s, image reconstruction from incomplete data remains challenging |[30l . Il33ll . 
Il42l . In this study, two data completion strategies were investigated Ii30l . [[33], l!42l to synthesize 
a complete data set, from which the WISE method could be effectively applied. 

One strategy was to fill the missing data with pressure corresponding to a homogeneous 
medium as 


r combHl 

lm Jm rec L+Z 


fe:," cpl ] 


J-- 1 (m rec )L+Z> 


[g m Jm rec L+Z, 

for m rec = 0,1, • • • , M — 1, where 6 ]p Mi 


p-incpl r- 
’ Sm 


if m rec e M^ od 

otherwise, 

^ reci , and g“ mbH € 


(27) 


pML 


, denote the 


computer-simulated (with a homogeneous medium), the measured incomplete, and the combined 


January 5, 2015 


DRAFT 








16 


complete data vectors at the m-th data acquisition, respectively. The mapping J m l : M^ od fa 
{0,1, • • • , N iec — 1} denotes the inverse operator of J m as 




rec M—N Iec 

m — m — 


if 


M-N 1 ' 


m rec — m + 


< m rec — m < 


M+N ri 


2 ' z z (28) 
if - M-N™ < m rec_ m< -M+N™ . 


This data completion strategy is based on the assumption that the back-scatter from breast tissue 
in an appropriately sound speed-matched water bath is weak. This assumption suggests that the 
missing measurements can be replaced by the corresponding pressure data that would have been 
produced in the absence of the object. 

The second, more crude, data completion strategy was to simply fill the missing data with 


zeros, i.e., 


r combO 
[om 


m Iec L+l ~ 


0, otherwise, 


(29) 


where g"’"' 1 ’ 0 denotes the data completed with the second strategy. 

5) Bent-ray image reconstruction: A bent-ray method was also employed to reconstruct im¬ 
ages. Details regarding the time-of-flight estimation and algorithm implementation are provided 
in Appendix O 


V. Computer-simulation results 
A. Images reconstructed from idealized data 

The images reconstructed from the noise-free, non-attenuated, data by use of the WISE method 
with 199 iterations and the sequential waveform inversion method with 43 iterations are shown 
in Fig. [5]-(a) and (b). As expected [f23ll . Il54ll . both images are more accurate and possess higher 
spatial resolution than the one reconstructed by use of the bent-ray reconstruction algorithm 
displayed in Fig. 0-(c). Profiles through the reconstructed images are displayed in Fig. [6l The 
images shown in Fig. (5]-(a) and -(b) possess similar accuracies as measured by their root-mean- 
square errors (RMSEs), namely, 1.08 x HU 3 for the former and 1.19 x 10 ~ 3 for the latter. 
The RMSE was computed as the Euclidean distance between the reconstructed image and the 
sound speed phantom vector c, averaged by the 256 x 256 pixels of the ROI sketched in Fig. Q] 
However, the reconstruction of Fig. [5]- (a) required only about 1.7% of the computational time 
required to reconstruct Fig.[5]-(b), namely, 1.4 hours for the former and 81.4 hours for the latter 
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respectively. This is because the WISE method required only 1018 wave equation solver runs 
which is significantly less than the 57088 wave equation solver runs required by the sequential 
waveform inversion method. With a similar number of wave equation solver runs, (e.g., 1024), 
one can complete only a single algorithm iteration by use of the sequential waveform inversion 
method. The corresponding image, shown in Fig. O-fd), lacks quantitative accuracy as well as 
qualitative value for identifying features. The results suggest that the WISE method maintains 
the advantages of the sequential waveform inversion method while significantly reducing the 
computational time. 

B. Convergence of the WISE method 

Images reconstructed from noise-free, non-attenuated, data by use of the WISE method contain 
radial streak artifacts when the algorithm iteration number is less than 100, as shown in Figs. [7J- 
(a-c). Profiles through these images are displayed in [8} The streaks artifacts are likely caused by 
crosstalk introduced during the source encoding procedure ll3Tl . [ |4T| . However, these artifacts 
are effectively mitigated after more iterations as demonstrated by the image reconstructed after 
the 199-th iteration in Fig. [5]-(a) and its profile in Fig. [6j The quantitative accuracy of the 
reconstructed images is improved with more iterations as shown in Fig. [8} 

Figure [U-(a) reveals that the WISE method requires a larger number of algorithm iterations 
than does the sequential waveform inversion method to achieve the same RMSE. The RMSE 
of the images reconstructed by use of the WISE method appears to oscillate around 1.0 x 10 -3 
after the first 100 iterations while the sequential waveform inversion method can achieve a lower 
RMSE. However, as shown previously in Fig. [5]-(a) and the corresponding profile in Fig. [6} 
after additional iterations the image reconstructed by use of the WISE method achieves a high 
accuracy. Moreover, to achieve the same accuracy as the sequential waveform inversion method, 
the WISE method requires a computation time that is reduced by approximately two-orders 
of magnitude, as suggested by Fig. HHb). We also plotted the cost function value against the 
number of iterations in Fig. IU-(c). Note that for the WISE method, the cost function value was 
approximated by the current realization of ^||g w — MH c s w || 2 . These plots suggest that, in this 
particular case, the WISE method appears to approximately converge after 200 iterations. For 
example, the images reconstructed after 199 (Fig. 0-(a)) and 250 (Fig.[7]-(d)) iterations are nearly 
identical. 
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C. Images reconstructed from non-attenuated data containing noise 

Images reconstructed by use of the WISE method with a quadratic penalty and the WISE 
method with a TV penalty from noisy, non-attenuated, data are presented in Fig. [TOj All images 
were obtained after 1024 algorithm iterations. The WISE method with a quadratic penalty 
effectively mitigates image noise as shown in Figs. ITUl-fa-c). at the expense of image resolution, 
as expected. Figure [[ORd) shows an image reconstructed by use of the WISE method with a 
TV penalty. The image appears to possess a similar resolution but a lower noise level than the 
image in Fig. ITOkb) that was reconstructed by use of the WISE method with a quadratic penalty. 
We also compared the convergence rates of the WISE method and the sequential waveform 
inversion methods when both utlize a TV penalty and the same regularization parameter. As 
shown in Fig. QT1 the convergence properties of the penalized methods follow similar trends 
as the un-penalized methods, which were discussed above and shown in Fig. [9] Even though it 
required a larger number of algorithm iterations, the WISE method reduced the computation time 
by approximately two-orders of magnitude as compared to the sequential waveform inversion 
method. 

D. Images reconstructed from acoustically attenuated data 

Our current implementation of the WISE method assumes an absorption-free acoustic medium. 
This assumption can be strongly violated in practice. In order to investigate the robustness of 
the the WISE method to model errors associated with ignoring medium acoustic absorption, 
we applied the algorithm to the acoustically attenuated data that were produced as described in 
Section IIV-C1 As shown in Fig. O when acoustic absorption is considered, the amplitude of 
the measured pressure is attenuated by approximately a factor of 2. The wavefront (See Fig. fl2l - 
(a)) remains very similar to that when medium absorption is ignored (See Fig. |4]-(a)). Medium 
absorption has the largest impact on the pressure data received by transducers located opposite 
the emitter as shown in Fig. Q2Kb). The shape of the pulse profile remains very similar as shown 
in Fig. [H-(c) and -(d), suggesting that waveform dispersion may be less critical than amplitude 
attenuation in image reconstruction for this phantom. 

Images reconstructed by use of the WISE method with a TV penalty from noise-free and 
noisy attenuated data are shown in Figs. lT3l-(a) and (b). Image profiles are shown in Fig. IT3T-(c). 
Although these images contain certain artifacts that were not produced in the idealized data 
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studies, most object structures remain readily identified. These results suggest that the WISE 
method with a TV penalty can tolerate data inconsistencies associated with neglecting acoustic 
attenuation in the imaging model, at least to a certain level with regards to feature detection 
tasks. 

E. Images reconstructed from idealized incomplete data 

The wavefront of the noise- and attenuation-free pressure wavefield when the object is absent 
(Fig.lT4l-('aV) appears to be very similar to that when the object is present (Fig. |4]-(a)). As expected, 
the largest differences are seen in the signals received by the transducers located opposite of 
the emitter, as shown in Fig. [T4l-(b'). As seen in Fig. |T4]-(c), the time traces received by the 
40-th transducer are nearly identical when object is present and absent. This is because the 
back-scattered wavefield is weak for breast imaging applications. These results establish the 
potential efficacy of the data completion strategy of filling the missing data with the pressure 
data corresponding to a water bath. 

The image reconstructed from the measurements completed with pressure data corresponding 
to a water bath is shown in Figfl5l-(al. As revealed by the profile in Figfl5l-(cl. this image is 
highly accurate. Alternatively, the image reconstructed from the the data completed with zeros 
contains strong artifacts as shown in Fig. |T3]-(b). These results suggest that the WISE method can 
be adapted to reconstruct images from incomplete data, which is particularly useful for emerging 
laser-induced USCT imaging systems lfT3il - ffl5l . 

VI. Experimental validation 

A. Data acquisition 

Experimental data recorded by use of the SoftVue USCT scanner |[55ll was utilized to further 
validate the WISE method. The scanner contained a ring-shaped array of radius 110 mm that 
was populated with 2048 transducer elements. Each element had a center frequency of 2.75 
MHz, a pitch of 0.34 mm, and was elevationally focused to isolate a 3 mm thick slice of the 
to-be-imaged object. The transducer array was mounted in a water tank and could be translated 
with a motorized gantry in the vertical direction. Readers are referred to ll55l for additional 
details regarding the system. 
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The breast phantom was built by Dr. Ernie Madsen from the University of Wisconsin and pro¬ 
vides tissue-equivalent scanning characteristics of highly scattering, predominantly parenchymal 
breast tissue. The phantom mimics the presence of benign and cancerous masses embedded in 
glandular tissue, including a subcutaneous fat layer. Figure |T6] displays a schematic of one slice 
through the phantom. The diameter of the inclusions is approximately 12 mm. Table HQ presents 
the known acoustic properties of the phantom. 

During data acquisition, the breast phantom was placed near the center of the ring-shaped 
transducer array so that the distance between the phantom and each transducer was approximately 
the same. While scanning each slice, every other transducer element sequentially emits fan beam 
ultrasound signals towards the opposite side of the ring. The forward scattered and backscattered 
ultrasound signals are subsequently recorded by the same transducer elements. The received 
waveform was sampled at a rate of 12 MHz. The 1024 data acquisitions required approximately 
20 seconds in total. A calibration data set was also acquired in which the phantom object was 
absent. 

B. Data pre-processing 

48 bad channels were manually identified by visual inspection. After discarding these, the 
data set contained M = 976 acquisitions. Each acquisition contained N rec = 976 time traces. 
Each time trace contained L = 2112 time samples. The 976 good channels were indexed from 
0 to 975. The corresponding data acquisitions were indexed in the same way. A Hann-window 
low-pass filter with a cutoff frequency of 4 MHz was applied to every time trace in both the 
calibration and the measurement data. This data filtering was implemented to mitigate numerical 
errors that could be introduced by our second-order wave equation solver. 

C. Estimation of excitation pulse 

The shape of the excitation pulse was estimated as the time trace of the calibration data (after 
pre-processing) received by the 488-th receiver at the 0-th data acquisition. Note that the 488-th 
receiver was approximated located on the axis of the 0-th emitter, thus the received pulse was 
minimally affected by the finite aperture size effect of the transducers. Because our calibration 
data and measurement data were acquired using different electronic amplifier gains, the amplitude 
of the excitation pulse was estimated from the measurement data. More specifically, we simulated 
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the O-th data acquisition using the second-order pseudospectral k-space method and compared 
the simulated time trace received by the 300-th receiver with the corresponding measured time 
trace (after pre-processing). The ratio between the maximum values of these two traces was used 
to scale the excitation pulse shape. We selected the 300-th receiver because it resided out of the 
fan-region indicated in Fig. [H its received signals were unlikely to be strongly affected by the 
presence of the object. The estimated excitation pulse and its amplitude spectrum are displayed 
in Fig. [T7] Note that the experimental excitation pulse contained higher frequency components 
than did the computer-simulated excitation pulse shown in Fig. [3} 

D. Synthesis of combined data 

As discussed in Section IIV-C41 signals received by receivers located near the emitter can be 
unreliable ll23l . Our experimental data, as shown in Fig. ITSl-fa). contained noise-like measure¬ 
ments for the receivers indexed from 0 to 200, and from 955 to 975, in the case where the 0-th 
transducer functioned as the emitter. Also, our point-like transducer assumption introduces larger 
model mismatches for the receivers located near the emitter. As shown in Figs. HEKc) and -(d), 
even though the simulated time trace received by the 300-th receiver matches accurately with 
the experimentally measured one, the simulated time trace received by the 200-th receiver is 
substantially different compared with the experimentally measured one. In order to minimize the 
effects of model mismatch, we replaced these unreliable measurements with computer-simulated 
water bath data, as described in Section IIV-CI We designated the time traces received by the 
512 receivers located on the opposite side of the emitter as the reliable measurements for each 
data acquisition. The 0-th data acquisition of the combined data is displayed in Fig. [T8T-(b). 

E. Estimation of initial guess 

The initial guess for the WISE method was obtained by use of the bent-ray reconstruction 
method described in Appendix O We first filtered each time trace of the raw data by a band¬ 
pass Butterworth filter (0.5MHz - 2.5MHz). Subsequently, we extracted the TOF by use of the 
thresholding method with a thresholding value of 20% of the peak value of each time trace. 
The bent-ray reconstruction algorithm was applied for image reconstruction with a measured 
background sound speed 1.513 mm//is. The resulting image is shown in Fig. [I9l-(a) and has a 
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pixel size of 1 mm. Finally, the image was smoothed by convolving it with a 2D Gaussian kernel 
with a standard deviation of 2 mm. 

F. Image reconstruction 

We applied the WISE method with a TV penalty to the combined data set. The second-order 
wave equation solver was employed with a calculation domain of dimensions 512.0 x 512.0 
mm 2 . The calculation domain was sampled on a 2560 x 2560 Cartesian grid with a grid spacing 
of 0.2 mm. On a platform consisting of dual quad-core CPUs with a 3.30 GHz clock speed, 64 
GB RAM, and a single NVIDIA Tesla K20 GPU, each numerical solver run, took 40 seconds 
to calculate the pressure data for 2112 time samples. Knowing the size of the phantom, we set 
the reconstruction region to be within a circle of diameter 128 mm, i.e., only the sound speed 
values of pixels within the circle were updated during the iterative image reconstruction. We 
swept the value of /3 TV over a wide range to investigate its impact on the reconstructed images. 

G. Images reconstructed from experimental data 

As shown in Fig. [[H the spatial resolution of the image reconstructed by use of the WISE 
method with a TV penalty is significantly higher than that reconstructed by use of the bent-ray 
model-based method. In particular, the structures labeled ‘A’ and ‘B’ possess clearly-defined 
boundaries. This observation is further confirmed by the profiles of the two images shown in 
Fig. [20l In addition, the structure labeled ‘C’ in Fig. [T9l-(b) is almost indistinguishable in the 
image reconstructed by use of the bent-ray model-based method (see Fig. EU-(a)). The improved 
spatial resolution is expected because the WISE method takes into account high-order acoustic 
diffraction, which is ignored by the bent-ray method ll23l . Though not shown here, for the 
bent-ray method, we investigated multiple time-of-flight pickers If25l and systematically tuned 
the regularization parameter. As such, it is likely that Fig. [T9Kaf represents a nearly optimal 
bent-ray image in terms of the resolution. This resolution also appears to be similar to previous 
experimental results reported in the literature lf26l . 

The convergence properties of the WISE method with a TV penalty with experimental data 
were consistent with those observed in the computer-simulation studies. Images reconstructed by 
use of 10, 50, and 300 algorithm iterations are displayed in Fig. ED The image reconstructed by 
use of 10 iterations contains radial streak artifacts that are similar in nature to those observed in 
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the computer-simulation studies. These artifacts were mitigated after more iterations. The image 
reconstructed after 300 iterations (Fig. [2lT-(d')l appears to be similar to that after 200 iterations 
(Fig. [j9]-(b)), suggesting that the WISE method with a TV penalty is close to convergence after 
about 200 iterations. The time required to complete 200 iterations was approximately 14 hours. 
The estimated time it would take for the sequential waveform inversion method to produce 
a comparable image is approximately one month, assuming the same number of iterations is 
required as in the computer-simulation studies (i.e., 40). 

Despite the nonlinearity of the WISE method, the impact of the TV penalty appears to be 
similar to that observed in other imaging applications |[52l . If56)l (see Fig. [22]). Though not shown 
here, the impact of the quadratic penalty is also similar. As expected, a larger value of j3 reduced 
the noise level at the expense of spatial image resolution. These results suggest a predictable 
impact of the penalties on the images reconstructed by use of the WISE method. 

VII. Summary 

It is known that waveform inversion-based reconstruction methods can produce sound speed 
images that possess improved spatial resolution properties over those produced by ray-based 
methods. However, waveform inversion methods are computationally demanding and have not 
been applied widely in USCT breast imaging. In this work, based on the time-domain wave 
equation and motivated by recent mathematical results in the geophysics literature, the WISE 
method was developed that circumvents the large computational burden of conventional wave¬ 
form inversion methods. This method encodes the measurement data using a random encoding 
vector and determines an estimate of the sound speed distribution by solving a stochastic opti¬ 
mization problem by use of a stochastic gradient descent algorithm. With our current GPU-based 
implementation, the computation time was reduced from weeks to hours. The WISE method was 
systematically investigated in computer-simulation and experimental studies involving a breast 
phantom. The results suggest that the method holds value for USCT breast imaging applications 
in a practical setting. 

Many opportunities remain to further improve the performance of the WISE method. As shown 
in Fig. Q33 images reconstructed by use of the WISE method can contain certain artifacts that 
are not present in the image reconstructed by use of the bent-ray method. An example of such 
an artifact is the dark horizontal streak below the structure C. Because of the nonlinearity of the 
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image reconstruction problem, it is challenging to determine whether these artifacts are caused 
by imaging model errors or by the optimization algorithm, which might have arrived at a local 
minimum of the cost function. A more accurate imaging model can be developed to account 
for out-of-plane scattering, the transducer finite aperture size effect, acoustic absorption, as well 
as other physical factors. Also, the stochastic gradient descent algorithm is one of the most 
basic stochastic optimization algorithms. Numerous emerging optimization algorithms can be 
employed lf33ll . Kill to improve the convergence rate. In addition, there remains a great need to 
compare the WISE method with other existing sound speed reconstruction algorithms ||T9ll , ll40ll . 

There remains a need to conduct additional investigations of the numerical properties of the 
WISE method. Currently, a systematic comparison of the statistical properties of the WISE and 
the sequential waveform inversion method is prohibited by the excessively long computation 
times required by the latter method. This comparison will be interesting when a more efficient 
wave equation solver is available. Given the fact that waveform inversion is nonlinear and 
sensitive to its initial guess, it becomes important to investigate how to obtain an accurate 
initial guess. We also observed that the performance of the WISE method is sensitive to how 
strong the medium heterogeneities are and the profile of the excitation pulse. An investigation 
of the impact of the excitation pulse the numerical properties of the image reconstruction may 
help optimize hardware design. In addition, quantifying the statistics of the reconstructed images 
will allow application of task-based measures of image quality to be applied to guide system 
optimization studies. 


Appendix A 

Continuous-to-Discrete USCT Imaging Model 

In practice, each data function g m (r,t) is spatially and temporally sampled to form a data 
vector g m E M ArrecL , where N Tec and L denote the number of receivers and the number of time 
samples, respectively. We will assume that N Tec and L do not vary with excitation pulse. Let 
[g m]n™ c L+i denotes the (n rec L + Z)-th element of g m . When the receivers are point-like, g m is 
defined as 

[gm]n*w = 9 m(r(m , n rec ), /A 1 ), (30) 

where the indices n rec and / specify the receiver location and temporal sample, respectively, and 
A 1 is the temporal sampling interval. The vector r(m, n rec ) E 1 l m denotes the location of the 
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n rec -th receiver at the m-th data acquisition. 

A C-D imaging model for USCT describes the mapping of c(r) to the data vector g m and 
can be expressed as 


[§m]n rec L+£ AA-rnU 


r=r(m,7i rec ),t=ZA t 


for 


ec =o,i,--- ,jv reo —1 

i=0,l,- - - ,L—1 


(31) 


Note that the acousto-electrical impulse response If57l of the receivers can be incorporated into 
the C-D imaging model by temporally convolving s m (r, t ) in Eqn. ([!) with the receivers’ acousto¬ 
electrical impulse response if we assume all receiving transducers share an identical acousto¬ 
electrical impulse response. 


Appendix B 

Frechet derivative of data fidelity term 


Consider the integrated squared-error data misfit function, ll22l . If23l 

1 m 3 f r^ 9 

^ C (c) = - ^ / dr dt\gjj,t) -g m (r,t)] 2 , (32) 

m =0 

where g m (r,t) and g m (r,t) denote the measured data function and the predicted data function 
computed by use of Eqn. © with the current estimate of c(r). 

Both the sequential and WISE reconstruction method described in Section Hill require knowl¬ 
edge of the Frechet derivatives of J rCC (c) and lZ cc (c) with respect to c, denoted by V c J r( C 
and V c 1Z CC , respectively. The calculation of V c 7 Z cc can be readily accomplished for quadratic 
smoothness penalties If52ll . l(58l . For the integrated squared error data misfit function given in 
Eqn. (l32l) . V c J 700 can be computed via an adjoint state method as lf28ll . Il59l . Il60l 

X f T Q 2 

Vc ‘ 7rCC = S / dt( lm{r,T -t)—p m (r,t), (33) 

c 3 (r) ? ^ 0 Jo d 2 t 

where q m (r, t ) e L 2 (M 3 x [0, oo)) is the solution to the adjoint wave equation. The adjoint wave 
equation is defined as 


V 2 g m (r,t) - 


i d 2 


9m(r,t) = —T m (r, t), 


(34) 


c^(r) d 2 t 

where r m (r, t) = g m { r, T — t) — g m ( r, T — t ). The adjoint wave equation is nearly identical in 
form to the wave equation in Eqn. ([[]) except for the different source term on the right-hand 
side, suggesting the same numerical approach can be employed to solve both equations. Since 
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one needs to solve Eqns. CD) and (l34l) M times in order to calculate V c J r( C , it is generally true 
that the sequential waveform inversion is computationally demanding even for a 2D geometry 

ED. 


Appendix C 

Bent-ray model-based sound speed reconstruction 

We developed an iterative image reconstruction algorithm based on a bent-ray imaging model. 
The bent-ray imaging model assumes that an acoustic pulse travels along a ray path that connects 
the emitter and the receiver and accounts for the refraction of rays, also known as ray-bending, 
through an acoustically inhomogeneous medium. For each pair of receiver and emitter, the travel 
time, as well as the ray path, is determined by the medium’s sound speed distribution. Given 
the travel times for a collection of emitter-and-receiver pairs distributed around the object, the 
medium sound speed distribution can be iteratively reconstructed. This bent-ray model-based 
sound speed reconstruction (BRSR) method has been employed in the USCT literature Il26ll . 

m, ED- 

In order to perform the BRSR, we extracted a TOF data vector from the measured pressure 
data. Denoting the TOF data vector by T e R MNiec , each element of T represented the TOF 
from each emitter-and-receiver pair. The extraction of the TOF was conducted in two steps. First, 
we estimated the difference between the TOF when the object was present and the TOF when 
the object was absent by use of a thresholding method ll25l . Il64l . In particular, 20% of the peak 
value of each time trace was employed as the thresholding value. Second, a TOF offset was 
added to the estimated difference TOF for each emitter-and-receiver pair to obtain the absolute 
TOF, where the TOF offset was calculated according to the scanning geometry and the known 
background SOS. 

Having the TOF vector T, we reconstructed the sound speed by solving the following opti¬ 
mization problem: 

s = argrnin || T — K s s || 2 +/31Z(s), (35) 

S 

where s denotes the slowness (the reciprocal of the SOS) vector, and K s denotes the system 
matrix that maps the slowness distribution to the TOF data. The superscript ‘s’ indicates the 
dependence of K s on the slowness map. At each iteration, using the current estimate of the SOS, 
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a ray-tracing method ||65M was employed to construct the system matrix K s . Explicitly storing 
the system matrix in the sparse representation, we utilized the limited BFGS method ff66l to 
solve the optimization problem given in Eqn. (1351) . The estimated slowness was then converted 
to the sound speed by taking the reciprocal of s element-wisely. We refer the readers to If26ll . 
If62l - ll64l . If67il for more details about the BRSR method. 
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Tables 

TABLE I: Parameters of the numerical breast phantom Il44l - lf46l 


Structure 

index 

Tissue type 

Sound speed 

[mm-gs _1 ] 

Slope of attenuation 
[dB ■ (MHz) ~ v -cm -1 ] 

0 

Adipose 

1.47 

0.60 

1 

Parenchyma 

1.51 

0.75 

2 

Benign tumor 

1.47 

0.60 

3 

Benign tumor 

1.47 

0.60 

4 

Cyst 

1.53 

0.00217 

5 

Malignant tumor 

1.565 

0.57 

6 

Malignant tumor 

1.565 

0.57 

7 

Malignant tumor 

1.57 

0.57 
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TABLE II: Parameters of the experimental breast phantom 


Material 

Sound speed 

[mm-/rs 1 ] 

Attenuation coefficient 

at 2.5 MHz [dB/cm] 

Fat 

1.467 

0.48 

Parenchymal tissue 

1.552 

0.89 

Cancer 

1.563 

1.20 

Fibroadenoma 

1.552 

0.52 

Gelatin cyst 

1.585 

0.16 
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Figures 



Fig. 1: Schematic of a USCT system with a circular transducer array whose elements are 
indexed from 0 to 255. It shows the first data acquisition, where element-0 (in red) is emitting 
an acoustic pulse, while all 256 elements are receiving signals. The region-of-interest (ROI) is 
shaded in gray, and the dashed square box represents the physical dimensions (128 x 128 mm 2 ) 
of all reconstructed images. 
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Fig. 2: (a) Sound speed map [mm-/is 1 ] and (b) acoustic attenuation slope map 

[dB-(MHz) _?/ -cm~ 1 ] of the numerical breast phantom. 
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Fig. 3: (a) Normalized temporal profile and (b) amplitude spectrum of the excitation pulse 

employed in the computer-simulation studies. The dashed line in (b) marks the center frequency 
of excitation pulse at 0.82 MHz. 
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Fig. 4: Computer-simulated (a) noise-free and (b) noisy data vectors at the 0-th data acquisition, 
(c) Profiles of the pressure received by the 128-th transducer. The grayscale window for (a) and 
(b) is [-45,0] dB. 


January 5, 2015 


DRAFT 














38 



(a) (b) 



(c) (d) 

Fig. 5: Images reconstructed by use of (a) the WISE method after the 199-th iteration (1,018 
runs of the wave equation solver), (b) the sequential waveform inversion algorithm after the 
43-rd iteration (57, 088 runs of the wave equation solver), (c) the bent-ray model-based sound 
speed reconstruction method, and (d) the sequential waveform inversion algorithm after the 1-st 
iteration (1,024 runs of the wave equation solver) from the noise-free non-attenuated data. The 
grayscale window is [1.46,1.58] mm//xs. 
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Fig. 6: Profiles at y — 6.5 mm of the images reconstructed by use of the bent-ray TOF image 
reconstruction method and the WISE method from the noise-free non-attenuated data. 
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(a) (b) 



(c) (d) 

Fig. 7: Images reconstructed by use of the WISE method after (a) the 20-th, (b) the 50-th, (c) the 
100-th, and (d) the 250-th iteration from the noise-free, non-attenuated data set. The grayscale 
window is [1.46,1.58] mm//is. 
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Fig. 8: Profiles of the images reconstructed by use of the WISE method from the noise-free 
non-attenuated data after different numbers of iterations. 
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Fig. 9: Plots of the root-mean-square errors (RMSEs) of the images reconstructed from the 
noise-free data versus (a) the number of iterations and (b) the number of wave equation solver 
runs, (c) Plots of the cost function value versus the number of iterations. 
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Fig. 10: Images reconstructed from non-attenuated data contaminated with Gaussian random 
noise. Images (a-c) were reconstructed by use of the WISE method with a quadratic penalty 
with /3 q = 1.0 x 10" 3 , 1.0 x 10” 2 , and 1.0 x 10 respectively. Image (d) was reconstructed by 
use of the WISE method with a TV penalty with /3 TV = 5.0 x 10 -4 . The insert in the up right 
comer of each image is the zoomed-in image of the dashed black box, which contains 35 x 35 
pixels (17.5 x 17.5 mm 2 ). The grayscale window is [1.46,1.58] mm/ps. 
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Fig. 11: Plots of the root-mean-square errors (RMSEs) of the images reconstructed from the 
noisy data versus (a) the number of iterations and (b) the number of wave equation solver runs, 
(c) Plots of the cost function value versus the number of iterations. Both the WISE and the 
sequential waveform inversion methods employed a TV penalty with /3 TV = 5.0 x 10” 4 . 
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Fig. 12: (a) Computer-simulated noise-free attenuated pressure of the 0-th data acquisition. 

(b) The difference between the attenuated pressure data and the non-attenuated pressure data. 

(c) The temporal profiles and (d) the amplitude spectra of the pressure received by the 128-th 
transducer. The grayscale window for (a) and (b) is [—45,0] dB. 
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(a) (b) 



Fig. 13: (a) Image reconstructed by use of the WISE method from the noise-free attenuated data, 
(b) Image reconstructed by use of the WISE method with a TV penalty with /3 TV = 5.0 x 10~ 4 , 
from the noisy attenuated data. The grayscale window is [1.46,1.58] mm/fis. (c) Profiles at 
y = 6.5 mm of the reconstructed images. 
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Fig. 14: (a) Computer-simulated noise-free non-attenuated pressure data when the object is 
absent, (b) The difference between the pressure data when object is present and the pressure 
data when the object is absent, (c) Profiles of the pressure received by the 40-th transducer. The 
grayscale window for (a) and (b) is [—45,0] dB. 
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(a) (b) 



Fig. 15: Images reconstructed by use of the WISE method from noise-free combined data that 
are completed (a) with computer-simulated pressure corresponding to a homogeneous medium 
and (b) with zeros. The grayscale window is [1.46,1.58] mm//is. (c) Profiles at y = 6.5 mm of 
the images reconstructed by use of the WISE method from the two combined data sets. 
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Fibroadenoma 



Fig. 16: Schematic of the breast phantom employed in the experimental study. 


January 5, 2015 


DRAFT 


50 
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(a) 


(b) 


Fig. 17: (a) Normalized temporal profile and (b) amplitude spectrum of the excitation pulse 
employed in the experimental studies. The dashed line in (b) marks the center frequency of 
excitation pulse at 2.09 MHz. 


January 5, 2015 


DRAFT 










51 



4 


cd 



C/) 

C/) 

CD 


Q- 


-4 



Fig. 18: Zeroth acquisition of (a) the experimentally-measured raw data and (b) the combined 
data, respectively, and time traces at the 0-th acquisition received by (c) the 300-th receiver, and 
(d) the 200-th receiver, respectively. The grayscale window for (a) and (b) is [—45,0] dB. 
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(a) (b) 

Fig. 19: Images reconstructed from the experimentally measured phantom data by use of (a) the 
bent-ray model-based sound speed reconstruction method and (b) the WISE method with a TV 
penalty with (/3 TV = 1.0 x 10 2 ) after the 200-th iteration. The grayscale window is [1.49,1.57] 
mm//is. 
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Fig. 20: Profiles at (a) x = —24.0 mm and (b) x = 10.0 mm of the reconstructed images by 
use of the bent-ray model-based sound speed reconstruction method (light solid) and the WISE 
method with a TV penalty with /3 TV = 1.0 x 10 2 (dark dashed) from experimentally measured 
data. 
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(a) (b) 



(c) (d) 

Fig. 21: (a) The initial guess of the sound speed map and the images reconstructed by use of the 
WISE method with a TV penalty with (/3 TV = 1.0 x 10 2 ) after (b) the 10-th, (b) the 50-th and 
(d) the 300-th iteration, from the experimentally measured phantom data. The grayscale window 
is [1.49,1.57] mm//is. 
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(a) (b) 

Fig. 22: Images reconstructed by use of the WISE method with a TV penalty with (a) /3 TV = 
5.0 x 10 1 , and (b) /3 TV = 5.0 x 10 2 , from the experimentally measured phantom data. The 
grayscale window is [1.49,1.57] mm//js. 
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